# Manuscript replication file
# Santiago Lopez-Cariboni, 2022. "Political Regimes and Informal Social Insurance", Comparative Political Studies.

# ---- wd ----
wd <- '~/Dropbox/Research/Papers/ElectricityTheft/PoliticalRegimesLosses/replication_materials' 
setwd(wd)

# ---- packages ----
paquetes <- c("data.table", "foreign", 'panelView', 'fect', "lubridate", 
    "ggplot2", "dplyr", "mFilter", "texreg", "AER", "xts", "interflex", "neverhpfilter")
suppressPackageStartupMessages(sapply(paquetes,require,character.only=T))

# Plot the effects (use Johannes Karreth's function: https://github.com/jkarreth/JKmisc/blob/master/ggintfun.R)
# devtools::source_url("https://raw.githubusercontent.com/jkarreth/JKmisc/master/ggintfun.R")
source('code/functions/ggintfun.R')


# ---- data preparation ----
# Source data management file
# source('code/data_replication.R')
dt <- fread('data/dt_replication.csv')

# ---- Table 1 ----
m4.dev.c <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton 
    # + l.democracy
    + l.capacity
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[LDC==0,], na.action=na.omit)
nc4.dev.c <- length(unique(m4.dev.c$model[['as.factor(iso3c)']]))

# subset to developing countries for now on.
dt <- dt[LDC==1, ]

m1 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton
    + l.democracy
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
nc1 <- length(unique(m1$model[['as.factor(iso3c)']]))

m1dem <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy ==1,], na.action=na.omit)
nc1dem <- length(unique(m1dem$model[['as.factor(iso3c)']]))



m1aut <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
nc1aut <- length(unique(m1aut$model[['as.factor(iso3c)']]))


m2 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
nc2 <- length(unique(m2$model[['as.factor(iso3c)']]))

m3 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.capacity
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
nc3 <- length(unique(m3$model[['as.factor(iso3c)']]))


m4 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + l.capacity
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt, na.action=na.omit)
nc4 <- length(unique(m4$model[['as.factor(iso3c)']]))


v.names <- c(
    "Losses Cycle$_{t-1}$",
    "GDP Output Gap$_{t-1}$",
    "Democracy$_{t-1}$",
    "Democracy$_{t-1}$ $\\times$ GDP Output Gap$_{t-1}$",
    "State Capacity$_{t-1}$",
    "Imports$_{t-1}$",
    "Exports$_{t-1}$",
    "Population (log)$_{t-1}$",
    "Real GDP per capita (log)$_{t-1}$",
    "Electricity Consumption$_{t-1}$",
    "Population Density$_{t-1}$"
    )


texreg(list(m1dem,m1aut,m1,m2,m3,m4,m4.dev.c),
    file= "tables/2wfe.tex",
    label="tab:twfe",
    custom.header = list("Developing countries" =1:6),
    custom.model.names = c('(Democracies)','(Autocracies)', '(All)',   '(All)', '(All)', '(All)', "Developed \n countries"),
    caption="Ciclicality of Electricity Losses in Autocracies and Democracies (1971 - 2014)",
    dcolumn = TRUE,
    no.margin=FALSE,
    fontsize="scriptsize",
    single.row=FALSE,
    use.packages=FALSE,
    booktabs = TRUE,
    digits=2,
    float.pos="htbp",
    sideways=FALSE,
    omit.coef= "(year)|(iso3c)|(Intercept)",
    include.rsquared = FALSE,
    stars = c(0.01, 0.05, 0.1),
    custom.gof.rows = list("Year dummies" = c("Yes", "Yes", "Yes","Yes","Yes", "Yes", "Yes"), "Country fixed-effects" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
    	"Number of countries"= c(nc1dem, nc1aut, nc1, nc2, nc3, nc4, nc4.dev.c))
    # ,
    # reorder.gof = c(1, 2, 5, 3, 4),
    # custom.coef.names=v.names
    )



# ---- Table 2 ----
# 2SLS Results: Insturmenting for Democracy and GDP Output Gap
iv_dem <- ivreg(outgap.tdl ~ l.outgap.tdl
    + l.outgap.gdp.hamilton * l.democracy
    + as.factor(iso3c)
    + as.factor(year)
    | .
    - l.outgap.gdp.hamilton*l.democracy
    + l.outgap.gdp.hamilton*l.wreg.democracy
    ,
    data=dt, na.action=na.omit)
nciv_dem <- round(length(unique(iv_dem$model[['as.factor(iso3c)']])),2)
Fdem1 <- round(summary(iv_dem,  diagnostics = TRUE)$diagnostics[1,3],2)
Fint1 <- round(summary(iv_dem,  diagnostics = TRUE)$diagnostics[2,3],2)

iv_out <- ivreg(outgap.tdl ~ l.outgap.tdl
    + l.outgap.gdp.hamilton * l.democracy
    + as.factor(iso3c)
    + as.factor(year)
    | .
    - l.outgap.gdp.hamilton*l.democracy
    + l.w.outgap.gdp.hamilton*l.democracy
    ,
    data=dt, na.action=na.omit)# summary(iv_out,  diagnostics = TRUE)
nciv_out <- round(length(unique(iv_out$model[['as.factor(iso3c)']])),2)
Fog2 <- round(summary(iv_out,  diagnostics = TRUE)$diagnostics[1,3],2)
Fint2 <- round(summary(iv_out,  diagnostics = TRUE)$diagnostics[2,3],2)


iv_both <- ivreg(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * l.democracy
    + as.factor(iso3c)
    + as.factor(year)
    | .
    - l.outgap.gdp.hamilton*l.democracy
    + l.w.outgap.gdp.hamilton*l.wreg.democracy
    ,
    data=dt, na.action=na.omit)
# summary(iv_both,  diagnostics = TRUE)
nciv_both <- length(unique(iv_both$model[['as.factor(iso3c)']]))
Fog3 <- round(summary(iv_both,  diagnostics = TRUE)$diagnostics[1,3],2)
Fdem3 <- round(summary(iv_both,  diagnostics = TRUE)$diagnostics[2,3],2)
Fint3 <- round(summary(iv_both,  diagnostics = TRUE)$diagnostics[3,3],2)

iv_both_controls <- ivreg(outgap.tdl ~ l.outgap.tdl
    + l.outgap.gdp.hamilton * l.democracy
    + l.capacity    
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
    + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
    | .
    - l.outgap.gdp.hamilton*l.democracy
    + l.w.outgap.gdp.hamilton*l.wreg.democracy
    ,
    data=dt, na.action=na.omit)

nciv_both_controls <- length(unique(iv_both_controls$model[['as.factor(iso3c)']]))
Fog4 <- round(summary(iv_both_controls,  diagnostics = TRUE)$diagnostics[1,3], 2)
Fdem4 <- round(summary(iv_both_controls,  diagnostics = TRUE)$diagnostics[2,3], 2)
Fint4 <- round(summary(iv_both_controls,  diagnostics = TRUE)$diagnostics[3,3], 2)


v.names <- c(
    "Losses Cycle$_{t-1}$",
    "GDP Output Gap$_{t-1}$",
    "Democracy$_{t-1}$",
    "Democracy$_{t-1}$ $\\times$ GDP Output Gap$_{t-1}$",
    "State Capacity$_{t-1}$",
    "Imports$_{t-1}$",
    "Exports$_{t-1}$",
    "Population (log)$_{t-1}$",
    "Real GDP per capita (log)$_{t-1}$",
    "Electricity Consumption$_{t-1}$",
    "Population Density$_{t-1}$"
    )


texreg(list(iv_dem,iv_out,iv_both,iv_both_controls),
    file= "tables/IV2wfe.tex",
    label="tab:IVtwfe",
    custom.model.names = c('Democracy IV','OutputGap IV', 'Both IV', 'Both IV'),
    caption="2SLS Results: Insturmenting for Democracy and GDP Output Gap",
    dcolumn = TRUE,
    no.margin=FALSE,
    fontsize="scriptsize",
    single.row=FALSE,
    use.packages=FALSE,
    booktabs = TRUE,
    digits=2,
    float.pos="htbp",
    sideways=FALSE,
    omit.coef= "(year)|(iso3c)|(Intercept)",
    include.rsquared = FALSE,
    stars = c(0.01, 0.05, 0.1),
    custom.gof.rows = list(
    	"Democracy Insturment F-stat.:" = c(Fdem1, "", Fdem3, Fdem4),
    	"Ouput Gap Insturment F-stat.:" = c("", Fog2, Fog3, Fog4),
    	"Ouput Gap X Democracy Insturment F-stat.:" = c(Fint1,Fint2,Fint3,Fint4),
    	"Year dummies" = c("Yes", "Yes","Yes","Yes"), 
    	"Country fixed-effects" = c("Yes", "Yes", "Yes", "Yes"),
    	"Number of countries"= c(nciv_dem,nciv_out,nciv_both,nciv_both_controls)),
    reorder.gof = c(1, 2, 3, 4, 5, 6, 8, 7)
   # custom.coef.names=v.names
    )





# ---- Table 3 ----
# Ciclicality of Electricity Losses in Autocracies Moderated by Regime Duration
dt$lndura <- log(1+ dt$gwf_duration)


m_dura <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * lndura
    # + l.outgap.gdp.hamilton * l.democracy
    # + l.capacity
    # + l.imports
    # + l.exports
    # + log(l.pop)
    # + log(l.gdp.pcap)
    # + l.elec.cons.pc
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
summary(m_dura,  diagnostics = TRUE)
ncdura <- length(unique(m_dura$model[['as.factor(iso3c)']]))



m_dura2 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * lndura
    + l.capacity
    # + l.imports
    # + l.exports
    # + log(l.pop)
    # + log(l.gdp.pcap)
    # + l.elec.cons.pc
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
summary(m_dura2,  diagnostics = TRUE)
ncdura2 <- length(unique(m_dura2$model[['as.factor(iso3c)']]))

m_dura3 <- lm(outgap.tdl ~ l.outgap.tdl 
    + l.outgap.gdp.hamilton * lndura
    + l.capacity
    + l.imports
    + l.exports
    + log(l.pop)
    + log(l.gdp.pcap)
    + l.elec.cons.pc
        + l.pop.density
    + as.factor(iso3c)
    + as.factor(year)
     ,
    data=dt[democracy==0,], na.action=na.omit)
summary(m_dura3,  diagnostics = TRUE)
ncdura3 <- length(unique(m_dura3$model[['as.factor(iso3c)']]))


v.names <- c(
    "Losses Cycle$_{t-1}$",
    "GDP Output Gap$_{t-1}$",
    "Regime Duration$_{t-1}$",
    "Duration$_{t-1}$ $\\times$ GDP Output Gap$_{t-1}$",
    "State Capacity$_{t-1}$",
    "Imports$_{t-1}$",
    "Exports$_{t-1}$",
    "Population (log)$_{t-1}$",
    "Real GDP per capita (log)$_{t-1}$",
    "Electricity Consumption$_{t-1}$",
    "Population Density$_{t-1}$"
    )


    # "Electricity Service Access$_{t-1}$",
    # "$/mathbf{W}_{t-1}^{imp.comp.}/mathbf{y}_{t-1}$",

texreg(list(m_dura, m_dura2, m_dura3),
    file= "tables/dura.tex",
    label="tab:dura",
    custom.model.names = c('(1)','(2)', '(3)'),
    caption="Ciclicality of Electricity Losses in Autocracies Moderated by Regime Duration",
    dcolumn = TRUE,
    no.margin=FALSE,
    fontsize="scriptsize",
    single.row=FALSE,
    use.packages=FALSE,
    booktabs = TRUE,
    digits=2,
    float.pos="htbp",
    sideways=FALSE,
    omit.coef= "(year)|(iso3c)|(Intercept)",
    include.rsquared = FALSE,
    stars = c(0.01, 0.05, 0.1),
    custom.gof.rows = list("Year dummies" = c("Yes", "Yes", "Yes"), 
        "Country fixed-effects" = c("Yes", "Yes", "Yes"),
        "Number of countries"= c(ncdura, ncdura2, ncdura3)),
    reorder.gof = c(1, 2, 5, 3, 4),
    custom.coef.names=v.names
    )



# ---- Figure 2 ----
pd <- rbind(dt[democracy ==0, .(outgap.tdl, Democracy = as.factor(0))], dt[democracy ==1, .(outgap.tdl, Democracy = 1)])
pd <- na.omit(pd) 
pdf("figures/Fig2_outcome.pdf", width = 8, height = 5)
ggplot(pd, aes(x = outgap.tdl, group = Democracy, fill=Democracy)) + 
geom_density(alpha = 0.2)+
 scale_fill_manual(values=c("red","blue")) +
 xlim(c(-1, 1)) +
 xlab("T&D Losses (deviation from the trend)") +
 theme_minimal()
dev.off()


# ---- Figure 3 ----

events <- xts(c("Dictatorship (BMR)", "Democracy (BMR)"), 
              as.Date(c("1973-01-01", "1985-01-01")))
addEventLines(events, srt=90, pos=2, on = 1)
pdf("figures/Fig3_outputgap_ury.pdf", width = 5, height = 5)
# par(mfrow = c(1, 1))
# One country
ctry.dat <- dt[iso3c == "URY", .(year = ymd(year, truncated = 2L), GDP_pc = 100*log(gdp.pcap))]
ctry.dat <- as.xts(ctry.dat)
gdp_filtered <- yth_filter(ctry.dat, h = 4, p = 2, output = c("x", "trend", "cycle"))
gdp_filtered$Output_Gap <- (gdp_filtered[,1]-gdp_filtered[,2])/gdp_filtered[,2]
plot(gdp_filtered[, 1:2], grid.col = "white", up.col = "darkgreen", dn.col = "darkred", legend.loc = "topleft",
    main = "URUGUAY", , yaxis.right = FALSE,
     panels = 'lines(gdp_filtered[,4], col="darkblue", type="h", on=NA, main= "GDP Output-Gap")')
dev.off()

ctry.dat <- dt[iso3c == "URY", .(year = ymd(year, truncated = 2), 
    Output_Gap =  outgap.gdp.hamilton, 
    TDLosses =  outgap.tdl,
    PolityIV = polity2,
    Democracy = democracy)]
ctry.dat <- as.xts(ctry.dat)
plot(ctry.dat[, 2], grid.col = "white", up.col = "darkgreen", dn.col = "darkred", 
    legend.loc = "bottomright",
    main = "", type = "l", yaxis.right = FALSE
 )
pdf("figures/Fig3_outputgap_ury2.pdf", width = 5, height = 5)
events <- xts(c("Dictatorship (BMR)", "Democracy (BMR)"), 
              as.Date(c("1973-01-01", "1985-01-01")))
addEventLines(events, srt=90, pos=2, on = 1)

dev.off()


pdf("figures/Fig3_outputgap_phl.pdf", width = 5, height = 5)
# par(mfrow = c(1, 1))
# One country
ctry.dat <- dt[iso3c == "PHL", .(year = ymd(year, truncated = 2L), GDP_pc = 100*log(gdp.pcap))]
ctry.dat <- as.xts(ctry.dat)
gdp_filtered <- yth_filter(ctry.dat, h = 4, p = 2, output = c("x", "trend", "cycle"))
gdp_filtered$Output_Gap <- (gdp_filtered[,1]-gdp_filtered[,2])/gdp_filtered[,2]
plot(gdp_filtered[, 1:2], grid.col = "white", up.col = "darkgreen", dn.col = "darkred", legend.loc = "topleft",
    main = "PHILIPPINES", , yaxis.right = FALSE,
     panels = 'lines(gdp_filtered[,4], col="darkblue", type="h", on=NA, main= "GDP Output-Gap")')
dev.off()

ctry.dat <- dt[iso3c == "PHL", .(year = ymd(year, truncated = 2), 
    Output_Gap =  outgap.gdp.hamilton, 
    TDLosses =  outgap.tdl,
    PolityIV = polity2,
    Democracy = democracy)]
ctry.dat <- as.xts(ctry.dat)
plot(ctry.dat[, 2], grid.col = "white", up.col = "darkgreen", dn.col = "darkred", 
    legend.loc = "bottomright",
    main = "", type = "l", yaxis.right = FALSE
 )
pdf("figures/Fig3_outputgap_phl2.pdf", width = 5, height = 5)
events <- xts(c("Democracy (BMR)"), 
              as.Date(c("1986-01-01")))
addEventLines(events, srt=90, pos=2)
dev.off()



# ---- Figure 4 ----
# Effects of the Output-Gap on Electricity Losses conditional on Regime Type
pdf("figures/Fig4_effects.pdf", width = 8, height = 6)
ggintfun(obj = m2, varnames = c("l.democracy", "l.outgap.gdp.hamilton"), 
         varlabs = c("Democracy", "GDP Output Gap"),
         title = FALSE, rug = TRUE,
         twoways = TRUE)
dev.off()


# ---- Figure 5 ----
# Effects of the Output-Gap on Electricity Losses conditional on Autocratic Regime Duration
pdf("figures/Fig5_dura.pdf", width = 8, height = 6)
ggintfun(obj = m_dura, varnames = c("lndura", "l.outgap.gdp.hamilton"), 
         varlabs = c("Regime Duration (log)", "GDP Output Gap"),
         title = FALSE, rug = TRUE,
         twoways = TRUE)
dev.off()







